{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/obss/BIOBSS/blob/main/examples/ecg_processing.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__BIOBSS - ECG Signal Processing__\n",
    "\n",
    "_This notebook includes guidelines to help using BIOBSS for ECG signal processing and feature extraction._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__To run this notebook in Colab:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "git clone https://github.com/obss/biobss.git\n",
    "cd biobss\n",
    "pip install ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -U matplotlib\n",
    "import shutil\n",
    "from biobss import FIXTURES_ROOT\n",
    "shutil.move(\"/content/biobss/sample_data\",FIXTURES_ROOT.parent)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Import required packages:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import biobss\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import neurokit2"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of Contents\n",
    "1. [ECG Sample Data](#sampledata)<br>\n",
    "2. [ECG Signal Preprocessing](#ecg_pre)<br>\n",
    "    2.1. [Filtering](#ecg_filter)<br>\n",
    "    2.2. [Peak Detection](#ecg_peak)<br>\n",
    "    2.3. [Delineation](#ecg_waves)<br>\n",
    "    2.4. [Plotting](#ecg_plot)<br>\n",
    "4. [ECG Signal Quality Assessment](#ecg_sqa)<br>\n",
    "    4.1. [Clipping and Flatline Detection](#ecg_cf)<br>\n",
    "    4.2. [Physiological Limits](#ecg_pm)<br>\n",
    "    4.3. [Template Matching](#ecg_tm)<br>\n",
    "    4.4. [Simultaneous SQA Assessment](#sqa_all)<br>\n",
    "5. [ECG Feature Extraction](#ecg_features)<br>\n",
    "    5.1. [Features from R peaks](#ecg_f_rpeaks)<br>\n",
    "    5.2. [Features from P, Q, R, S, T waves](#ecg_f_waves)<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### __ECG Sample Data__\n",
    "<a id=\"sampledata\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ECG sample data is provided as a csv file in BIOBSS\\sample data. The data file contains an ECG signal of 5 minutes length, sampled at 256 Hz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load the sample data\n",
    "data, info = biobss.utils.load_sample_data(data_type='ECG')\n",
    "sig = np.asarray(data['ECG'])\n",
    "fs = info['sampling_rate']\n",
    "L = info['signal_length']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### __ECG Signal Preprocessing__\n",
    "<a id=\"ecg_pre\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### __Filtering__\n",
    "<a id=\"ecg_filter\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BIOBSS provides a filtering function which uses Scipy. This function can be used to implement Butterworth filter by defining the filter parameters (filter type, filter order, cutoff frequencies) as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Filter ECG signal by defining the filter parameters\n",
    "f_sig= biobss.preprocess.filter_signal(sig=sig, sampling_rate=fs, filter_type='bandpass',N=2,f_lower=0.5,f_upper=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an alternative, predefined filters can be used for each signal type. In order to use this option for ECG signal, signal_type should be selected as 'ECG'. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Filter ECG signal by using predefined filters\n",
    "filtered_ecg=biobss.preprocess.filter_signal(sig, sampling_rate=fs, signal_type='ECG', method='pantompkins')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### __Peak Detection__\n",
    "<a id=\"ecg_peak\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BIOBSS provides a specialized peak detection function for ECG signal. The function uses __ecgdetectors__ package and returns a dictionary of R-peak locations and amplitudes. The available methods are: 'pantompkins', 'hamilton' and 'elgendi'. \n",
    "For more, see: https://github.com/berndporr/py-ecg-detectors/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Detect peaks using 'pantompkins' method.\n",
    "locs_peaks=biobss.ecgtools.ecg_detectpeaks(sig,fs,'pantompkins')\n",
    "peaks = sig[locs_peaks]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### __Delineation__\n",
    "<a id=\"ecg_waves\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since ECG signal has a complex waveform including P wave, QRS complex and T wave; locations of other fiducial points may be required for calculation of some features. BIOBSS does not provide an ECG delineation function so __Neurokit2__ package can be used for this purpose. Neurokit's ___ecg_delineate___ function takes the signal and R-peak locations as input and returns a dictionary of fiducial locations. The function has three options for the delineation method, which are \"peak\" for a peak-based method, \"cwt\" for continuous wavelet transform or \"dwt\" (default) for discrete wavelet transform. The returned fiducials differ depending on the method selected for delineation. In this example, the method 'neurokit2' is selected and the fiducial locations can be accessed using the keys given below.\n",
    "\n",
    "__Reference__: Makowski, D., Pham, T., Lau, Z. J., Brammer, J. C., Lespinasse, F., Pham, H., Schölzel, C., & Chen, S. H. A. (2021). NeuroKit2: A Python toolbox for neurophysiological signal processing. Behavior Research Methods, 53(4), 1689–1696. https://doi.org/10.3758/s13428-020-01516-y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Delineate ECG signal using 'neurokit2' package. \n",
    "\n",
    "_, fiducials = neurokit2.ecg_delineate(ecg_cleaned=sig, rpeaks=locs_peaks, sampling_rate=fs, method='peak')\n",
    "\n",
    "p_peaks_locs = fiducials['ECG_P_Peaks']\n",
    "q_peaks_locs = fiducials['ECG_Q_Peaks']\n",
    "s_peaks_locs = fiducials['ECG_S_Peaks']\n",
    "t_peaks_locs = fiducials['ECG_T_Peaks']\n",
    "p_onset_locs = fiducials['ECG_P_Onsets']\n",
    "t_offset_locs = fiducials['ECG_T_Offsets']"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### __Plotting__\n",
    "<a id=\"ecg_plot\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "BIOBSS provides plotting functions specific to each signal type. In order to plot ECG signals, ___plot_ecg___ function can be used. The plots can be generated either using __Matplotlib__ or __Plotly__. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The signals and peaks should be provided as dictionaries, and the keys should be selected as shown below. To plot peaks, the parameter __show_peaks__ should be True and the parent keys of peaks dictionary should match with the keys of signals dictionary. In order to plot peaks for only some of the signals, it is enough to skip peak information in the peaks dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Generate inputs as dictionaries. The peaks are plotted for raw signal and skipped for filtered signal.\n",
    "signals={'Raw': sig, 'Filtered': filtered_ecg}\n",
    "peaks={'Raw': {'Ppeaks':p_peaks_locs, 'Qpeaks': q_peaks_locs, 'Rpeaks': locs_peaks}}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot ECG signal\n",
    "biobss.ecgtools.plot_ecg(signals=signals, peaks=peaks, show_peaks=True, figsize=(11, 7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot ECG signal using Plotly\n",
    "biobss.ecgtools.plot_ecg(signals=signals, peaks=peaks, method='plotly', show_peaks=True, width=800, height=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### __ECG Signal Quality Assessment__\n",
    "<a id=\"ecg_sqa\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Signal quality assessment is an important step which may be required prior to signal analysis. For this purpose, rule-based or machine learning based approaches can be used based on the morphological information from ECG waveform. BIOBSS has modules for assessing signal quality, which can be applied seperately or consecutively."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A clipped and flatline-including version of the ECG sample data can be generated to be used in the following steps. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Generate a modified sample data\n",
    "\n",
    "sig_modified = sig[:640].copy()\n",
    "#Clip the signal\n",
    "sig_modified[sig_modified > 0.2]=0.2\n",
    "sig_modified[sig_modified < -0.2]=-0.2\n",
    "#Add flatlines\n",
    "sig_modified[53:69] =sig_modified[53]\n",
    "sig_modified[217:233] =sig_modified[217]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the modified signal\n",
    "signals={'Raw': sig_modified}\n",
    "biobss.ecgtools.plot_ecg(signals, method='plotly', show_peaks=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### __Clipping and Flatline Detection__\n",
    "<a id=\"ecg_cf\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ___detect_flatline_clipping___ function from __sqatools__ subpackage is used to detect clipped or flat segments of a signal by setting the parameters and returns a dictionary of boundaries. \n",
    "\n",
    "Clipping occurs because of exceeding voltage limits of the signal conditioning circuits. Clipped segments can be detected by setting a clipping threshold and applying rules on the signal amplitude. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Detect clipped segments of ECG signal by setting a threshold value for signal amplitude\n",
    "clipped_segments=biobss.sqatools.detect_clipped_segments(sig_modified,threshold_pos=0.2,threshold_neg=-0.2)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Flatline detection differs from clipping such that it can occur at any level of signal amplitude and duration of flat segments also matters to detect flatlines. Thus, both an amplitude and duration threshold should be defined to apply the rules. Note that, amplitude threshold is set for the the minimum level of amplitude change required for a signal to be considered as non-flat segment, rather than absolute signal amplitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Detect flat segments of ECG signal by setting thresholds for signal amplitude change and duration of flat segments\n",
    "flatline_segments=biobss.sqatools.detect_flatline_segments(sig_modified,change_threshold=0.001, min_duration=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### __Physiological Limits__\n",
    "<a id=\"ecg_pm\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ___check_phys___ function can be used to check for physiological limits. This function takes R-peak locations as argument and returns a dictionary of boolean results of the applied rules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Peak detection\n",
    "locs_peaks=biobss.ecgtools.ecg_detectpeaks(sig_modified,fs,'pantompkins')\n",
    "peaks = sig_modified[locs_peaks]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the phsyiological limits can be compared to the accepted values in the literature. \n",
    "Reference: https://link.springer.com/book/10.1007/978-3-319-68415-4\n",
    "\n",
    "These values are defined as constants in the corresponding modules as given below. \n",
    "\n",
    "Physiological limits:\n",
    "- HR_MIN = 40 #minimum heart rate\n",
    "- HR_MAX = 180 #maximum heart rate\n",
    "- PP_MAX = 3 #maximum peak to peak interval\n",
    "- MAX_PP_RATIO = 2.2 #maximum P-P interval / minimum P-P interval \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check for physiological limits\n",
    "info=biobss.sqatools.check_phys(locs_peaks,fs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### __Template Matching__\n",
    "<a id=\"ecg_tm\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A common method for signal quality assessment is \"Template Matching\". This method is based on the expectation that pulse waveforms in an ECG segment will have similar morphology. A template is generated by aligning the pulses by their R-peaks and averaging them. Then, similarity of each pulse with the template is calculated using a measure. The ___template_matching___ function uses Pearson correlation as similarity measure. A threshold should be set for correlation coefficient below which the pulse is \"unacceptable\". The default value is set as 0.9 for the threshold. The function returns a tuple of correlation coefficients and a boolean result for the quality of the segment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Template matching \n",
    "info=biobss.sqatools.template_matching(sig_modified,locs_peaks)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### __Simultaneous SQA Assessment__\n",
    "<a id=\"sqa_all\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The methods defined above can also be applied simultaneously using the function ___sqa_ecg___ from __ecgtools__ subpackage. Methods can be provided as a list and the function returns a dictionary of results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info=biobss.ecgtools.sqa_ecg(ecg_sig=sig_modified, sampling_rate=fs, methods=['clipping', 'flatline','physiological','template'], threshold_pos=0.2, change_threshold=0.01, min_duration=10, peaks_locs=locs_peaks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### __ECG Feature Extraction__\n",
    "<a id=\"ecg_features\"></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ECG signal is mostly used for heart rate variability analysis and arrhythmia classification in the literature. For arrhythmia classification, morphological/time domain features are commonly used to train the machine learning models. BIOBSS package provides a set of functions to calculate some of the common ECG features in the literature. Some features can be calculated using only R-peak locations while others requires information on the other fiducial points (P, Q, S and T waves). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### __Features from R peaks__\n",
    "<a id=\"ecg_f_rpeaks\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ___from_Rpeaks___ function calculates morphological features using R-peak locations. These features are:\n",
    "- 'a_R': Amplitude of R peak\n",
    "- 'RR0': Previous RR interval\n",
    "- 'RR1': Current RR interval\n",
    "- 'RR2': Subsequent RR interval\n",
    "- 'RRm': Mean of RR0, RR1 and RR2\n",
    "- 'RR_0_1': Ratio of RR0 to RR1\n",
    "- 'RR_2_1': Ratio of RR2 to RR1\n",
    "- 'RR_m_1': Ratio of RRm to RR1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Calculate features from R peaks\n",
    "features_rpeaks = biobss.ecgtools.ecg_features.from_Rpeaks(sig, locs_peaks, fs)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to return segment-based (averaged) features, the parameter 'average' should be set as True."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Calculate features from R peaks\n",
    "features_rpeaks = biobss.ecgtools.ecg_features.from_Rpeaks(sig, locs_peaks, fs, average=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### __Features from P, Q, R, S, T waves__\n",
    "<a id=\"ecg_f_waves\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ___from_waves___ function calculates morphological features using locations of P, Q, R, S, T waves. These features are:\n",
    "- 't_PR': Time between P and R peak locations\n",
    "- 't_QR': Time between Q and R peak locations\n",
    "- 't_SR': Time between S and R peak locations\n",
    "- 't_TR': Time between T and R peak locations\n",
    "- 't_PQ': Time between P and Q peak locations\n",
    "- 't_PS': Time between P and S peak locations\n",
    "- 't_PT': Time between P and T peak locations\n",
    "- 't_QS': Time between Q and S peak locations\n",
    "- 't_QT':Time between Q and T peak locations\n",
    "- 't_ST': Time between S and T peak locations\n",
    "- 't_PT_QS': Ratio of t_PT to t_QS\n",
    "- 't_QT_QS': Ratio of t_QT to t_QS\n",
    "- 'a_PQ': Difference of P wave and Q wave amplitudes\n",
    "- 'a_QR': Difference of Q wave and R wave amplitudes\n",
    "- 'a_RS': Difference of R wave and S wave amplitudes\n",
    "- 'a_ST': Difference of S wave and T wave amplitudes\n",
    "- 'a_PS': Difference of P wave and S wave amplitudes\n",
    "- 'a_PT': Difference of P wave and T wave amplitudes\n",
    "- 'a_QS': Difference of Q wave and S wave amplitudes\n",
    "- 'a_QT': Difference of Q wave and T wave amplitudes\n",
    "- 'a_ST_QS': Ratio of a_ST to a_QS\n",
    "- 'a_RS_QR': Ratio of a_RS to a_QR\n",
    "- 'a_PQ_QS': Ratio of a_PQ to a_QS\n",
    "- 'a_PQ_QT': Ratio of a_PQ to a_QT\n",
    "- 'a_PQ_PS': Ratio of a_PQ to a_PS\n",
    "- 'a_PQ_QR': Ratio of a_PQ to a_QR\n",
    "- 'a_PQ_RS': Ratio of a_PQ to a_RS\n",
    "- 'a_RS_QS': Ratio of a_RS to a_QS\n",
    "- 'a_RS_QT': Ratio of a_RS to a_QT\n",
    "- 'a_ST_PQ': Ratio of a_ST to a_PQ\n",
    "- 'a_ST_QT': Ratio of a_ST to a_QT"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If location information on all fiducials are not available or required, these can be skipped in the fiducials dictionary. Then, only features related to the given fiducials will be calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Calculate features from P, Q, R, S, T waves\n",
    "features_waves = biobss.ecgtools.ecg_features.from_waves(sig, locs_peaks, fiducials, fs)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to calculate segment-based (averaged) features, the parameter 'average' should be set as True."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Calculate features from P, Q, R, S, T waves\n",
    "features_waves = biobss.ecgtools.ecg_features.from_waves(sig, locs_peaks, fiducials, fs, average=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
